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We discuss the analytic computation of autocorrelation functions for the generalized Hybrid Monte Carlo algorithm 
applied to free field theory and compare the results with numerical results for the 0(4) spin model in two 
dimensions. We explain how the dynamical critical exponent z for some operators may be reduced from two to 
one by tuning the amount of randomness introduced by the updating procedure, and why critical slowing down 
is not a problem for other operators. 



1. GENERALIZED HMC 

The work reported here extends the results first 
presented in [1]. We begin by recalling that a 
Markov Process will converge to some distribu- 
tion of configurations if it is constructed out of 
update steps each of which has the desired distri- 
bution as a fixed point, and which taken together 
are ergodic. The generalized HMC algorithm is 
constructed out of two such steps. 

1.1. Molecular Dynamics Monte Carlo 

This consists of three parts: (1) MD: an ap- 
proximate integration of Hamilton's equations on 
phase space which is exactly area-preserving and 
reversible; U(t) : (</>,tt) h- y tt') where det U = 
1 and U(t) = U{— t) -1 . (2) A momentum flip 
F : ir I— > —ir. (3) MC: a Metropolis accept/reject 
test. 

1.2. Partial Momentum Refreshment 

This mixes the Gaussian-distributed momenta ir 
with Gaussian noise £: 

/ tt' \ _ / cos 9 sin 9 \ / 7T \ 
^ (' J ~ \ - sin 9 cos 9 J \ C J 

The HMC algorithm is the special case where 
9 = -|. 9 = corresponds to an exact version 
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of the MD or microcanonical algorithm (which is 
in general non-ergodic). The L2MC algorithm of 
Horowitz [2,3] corresponds to choosing arbitrary 
9 but MDMC trajectories of a single leapfrog in- 
tegration step. 

2. LEAPFROG EVOLUTION 

Consider a system of harmonic oscillators {<f> p } 
for p G TLy . The Hamiltonian on phase space is 

H = -| X^peZv i^p w p^ 2 )- ^ s ^ e Hamiltonian 
is diagonal we shall temporarily suppress the in- 
dex p and set u> p = 1. 

2.1. Evolution operator 

The leapfrog discretization of Hamilton's equa- 
tions can be written as 

The most general area-preserving reversible linear 
mapping may be parameterized as 

( cos[«(*r)*r] ™ f ;<ff r1 \ 

\ — p(Sr) sin[/«(Jr)Jr] cos[k(St)St] ) 

where k and p are even functions of St. For 
leapfrog we find that k = 1 + St 2 + St 4 + 
0(St 6 ) and p = 1 - ISt 2 - ^8t a + 0(St 6 ). 
With this parameterization it is easy to see that 
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the evolution operator for a trajectory of t j St 
leapfrog steps is 

\k(Si 



2D 0(4) model, HMC, x=1 .0, 64 Lattice p=2.0 



u(t) = I cos K^k] ^ -^j 

\ — p(Sr) sin [«(Jr)r] cos[k(St)t] 



3. ACCEPTANCE RATES 

Given the evolution matrix we can compute the 
average Metropolis acceptance rate (P* cc ). The 
probability distribution of SH may be evaluated 
using Laplace's method to give an asymptotic ex- 
pansion in the lattice volume V 
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exp 



4<J£T> 



just as we would expect from the central limit 
theorem. The average Metropolis acceptance rate 
is thus 



(P 



.„) ~ erfc = erfc (^f^ 



The acceptance rate is a function of the "scaling" 
VSt 4 : The explicit dependence is 



variable x 



x v-> 



PGZ, 



(sin w p r) 2 w^. 



4. AUTOCORRELATION FUNCTIONS 



Markov Processes 

(/>jv) be a sequence of field con- 



4.1. 

Let {<j)i,<j)2, 

figurations generated by an equilibrated ergodic 
Markov process, and let (£!((/>)} denote the expec- 
tation value of some connected operator SI. We 
may define an unbiased estimator £1 over the finite 
sequence of configurations by SI = ^2t=i ^(<^t); 



As usual, we define C'n(£) 

autocorrelation function for SI 
the estimator SI is 



as the 
The variance of 



1 + 



I 



N 



where Aq = X^fci On(^) is the integrated auto- 
correlation function for the operator SI and N exp 
is the exponential autocorrelation time. This re- 
sult tells us that on average 1 + 2Aq correlated 
measurements are needed to reduce the variance 
by the same amount as a single truly independent 
measurement. 
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Figure 1. Comparison of acceptance rates for 2D 
0(4) model with 2D free field theory; n is the 
order of the integration scheme, with n = for 
leapfrog. The higher order integration schemes 
have a scaling variable of the form x = VSt 4 " +4: . 



4.2. Cost 

We shall calculate Aq as a function of MD time. 
To a good approximation the cost (£ of the com- 
putation is proportional to the total MD time for 
which we have to integrate Hamilton's equations. 
The cost per independent configuration is then 
<£ oc + 2 An)- The optimal trajectory length 
is obtained by minimizing the cost as a function 
of the parameters St, t, and 6 of the algorithm. 

4.3. Autocorrelation functions for polyno- 
mial operators 

We need to make a few simplifying assumptions: 

(1) The acceptance probability _P acc for each tra- 
jectory may be replaced by its average value; we 
neglect correlations in the acceptance probability 
between successive trajectories. Including such 
correlations leads to seemingly intractable com- 
plications. It is not obvious that our assumption 
corresponds to any systematic approximation ex- 
cept, of course, that it is valid when _P acc = 1. 

(2) (-P acc ) is assumed to be independent of trajec- 
tory length. This assumption is made purely for 
simplicity, as otherwise our results are expressed 
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in terms of particularly disgusting integrals. 
We may ignore the corrections of non-leading or- 
der in St to the MD evolution operator because 
for any given value of (-P acc ) there is a correspond- 
ing value of x which is 0(1), and thus St is of or- 
der V -1 / 4 . These corrections therefore only con- 
tribute to the autocorrelations through the accep- 
tance rate itself. 

We shall choose each trajectory r length inde- 
pendently from some distribution Pr(t), as this 
avoids the lack of ergodicity caused by choosing a 
fixed trajectory length which is a rational multi- 
ple of the period of any mode of the system. This 
is a disease of free field theory which in interact- 
ing models is removed to some extent by mode 
coupling. 

The integrated autocorrelation function for the 
connected squared magnetization M 2 with expo- 
nentially distributed trajectory lengths, Pn(t) = 
e-'/ r , is 

-4_P 2 cc w 2 r 2 cos 3 9 + 8_P acc w 2 r 2 cos 3 9- 
-4w 2 r 2 cos 3 9 + 2_P acc cos 3 9 - cos 3 9- 
-4P, cc uj 2 t 2 cos 2 9 - 4w 2 r 2 cos 2 9 - cos 2 9- 
-2P 2 cc uj 2 t 2 cos 9 - 6P, cc lo 2 t 2 cos 9+ 
+4uj 2 t 2 cos 9 - 2_P acc cos 9 - cos 9- 

-2P, cc i0 2 T 2 - 4co 2 t 2 - 1 

2 J P acc w 2 r 2 (cos 2 0-1) (P^ c cos 9 - cos 9 - 1) 

The integrated autocorrelation function for M 2 
with fixed length trajectories, Pn(t) = S(t — T), is 

-2_P acc cos 9 sin 2 9 + cos 9 sin 2 9- 
-2_P acc sin 2 {lot) sin 2 9 + sin 2 9+ 
+_P acc sin 2 (wr) cos 9 + _P acc sin 2 (wr) 

_P acc sin 2 (wr) (cos 9+1) sin 2 9 

5. RESULTS FOR C M ^(T) AND A M 2 

In order to understand the results let us consider 
the special case where 9 = (HMC) and _P acc = 1. 
In this case the Laplace transform of the autocor- 
relation function is 

r 3 + 2r 2 f3 + rj3 2 + 2m 2 r 
M 2 (P) - ^ + 2r/ j 2 + ^ 2 + 4m2 ^ + 2m2r , 

where we have written r = 1/r and u> = m. 



5.1. Autocorrelation functions 

The fact that Fm^{P) is a rational function in [3 
tells us that the autocorrelation function is a sum 
of exponentials. The exponents are the roots of 
the cubic denominator, and they are all real or 
one is real and the other two are complex conju- 
gates depending on the value of the mean trajec- 
tory length 1/r. 

5.2. Integrated autocorrelation function 

This is just A M 2 = F M 2{0) = 1 + i(^) 2 . 
Minimizing the cost by choosing r opt such that 
4^1 = we obtain r„„ t = l/v3m and 

thence Ajvfa (T opt ) = 5/2. The optimum trajectory 
length depends upon the operator being consid- 
ered. 

5.3. Dynamical critical exponent 

One of the most relevant measures of the effec- 
tiveness of an algorithm for studying continuum 
physics is the exponent z relating the cost (£ to 
the correlation length £ of the system. For free 
field theory the correlation length is just the in- 
verse mass, and thus we have z = 2 if we hold r 
constant, but z = 1 if we choose r = r opt which 
grows as l/m. 

5.4. Optimal choice for 9 

For the generalized HMC algorithm we should 
minimize the cost by varying both r and 9. The 
optimal choice of parameters when _P acc = 1 is 
to take 9^-0 and r — > 0, but this ignores the 
fact that the cost does not decrease when we take 
r smaller than the St required to obtain a rea- 
sonable Metropolis acceptance rate. If we choose 
r opt = St (L2MC) and the corresponding value 
for 9 opt we find that the cost is less than for the 
HMC case, but only by a constant factor. As the 
cost is only defined up to an implementation de- 
pendent constant factor anyhow we may conclude 
that generalized HMC does not appear to promise 
great improvements over HMC. 

5.5. L2MC 

Horowitz suggested that the L2MC algorithm was 
better than HMC not for the reasons discussed 
above, but because the acceptance rate is much 
higher for fixed St. Unfortunately, when 9 ^ -| 
we must flip the momenta upon rejecting a tra- 
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2D 0(4) model, 64 Lattice, p=2.0 
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Figure 2. Comparison of the HMC and L2MC 
algorithms. 



jectory, and this causes the system to perform a 
random walk unless the acceptance rate is very 
close to unity. Our data shows that in practice 
the cheapest solution appears to be HMC. 

6. RESULTS FOR A E 

So far all our results have been for the operator 
M 2 , but this is a rather special case as it couples 
solely to the slowest mode of the system. Let 
us now investigate the properties of the energy 
operator £=iE,^p- 

The Laplace transform of the connected auto- 
correlation function for the energy is Fe(/3) = 

i E , . The integrated 

autocorrelation function for the energy is 



At 



F E (0) = 1 + 



— V— = 1 + J- (- 1 

2t 2 ^w 2 2r 2,7m2 



and the optimal trajectory length is r OJ 



cr m2 1 ''/3) leading to an integrated autocorrela- 
tion function value of A£;(r opt ) = 5/2. 
In order to determine the dynamical critical ex- 
ponent z we need to evaluate the spectral sum 



6.1. Spectral sums for 2D free field theory 

Using Poisson resummation we find that in the 
thermodynamic limit 



V ^ 



v ^ - p 

Px.Py 

where a 



2\/ab 
tt(1 - ab)~ 



-K 



±m 2 + 3 



plete elliptic integral. For small m we find that 
a^~P « T-ln-^i-, and hence z = 0. Note that 
this does not mean that the cost does not increase 
with increasing £, but that it increases only loga- 
rithmically. 

7. CONCLUSIONS 

(1) Too little noise increases critical slowing down 
because the system is too weakly ergodic. (2) Too 
much noise increases critical slowing down be- 
cause the system takes a drunkard's walk through 
phase space. (3) To attain z = 1 for all operators 
(and especially for the exponential autocorrela- 
tion time) one must be able to tune the amount of 
noise suitably. (4) For operators such as E, crit- 
ical slowing down is unimportant because they 
only couple weakly to the slowest modes of the 
system. 
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V(im 2 + 3) 2 -l, b = 
and K(k) is a com- 



